home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / SoundAndMusic / cmix / lpc / analysis / alpole.c next >
C/C++ Source or Header  |  1991-12-10  |  3KB  |  146 lines
  1. #include "../../H/sfheader.h"
  2. #include <stdio.h>
  3. #include <sys/file.h>
  4. #include <sys/types.h>
  5. #include <sys/stat.h>
  6. #include <signal.h>
  7. #include <errno.h>
  8. #include <math.h>
  9.  
  10.  
  11. int NPOLE;
  12. int FRAME;
  13. int NPP1;
  14.  
  15. #define POLEMAX  100
  16. #define FLOAT 4
  17. #define FRAMAX   500 
  18. #define NDATA 4 /* number of data values stored with frame */
  19.  
  20. alpol(sig, errn, rms1, rms2, c)
  21. double *errn, *rms1, *rms2, *c;
  22. int *sig;
  23. {
  24.     double a[POLEMAX][POLEMAX], v[POLEMAX], b[POLEMAX];
  25.     double x[FRAMAX], y[FRAMAX];
  26.     double *vp=v, *bp=b, *xp=x, *yp=y;
  27.     double sum, sumx, sumy;
  28.     int k1, i, l, k, limit, j;
  29.  
  30.     for (xp=x; xp-x < FRAME ;++xp,++sig) 
  31.         *xp = *sig;
  32.     k1 = NPOLE + 1;
  33.     for (i=0; i < NPOLE ;++i)  {
  34.         sum = (double) 0.0;
  35.         for (k=NPOLE; k < FRAME ;++k)
  36.             sum += x[k-(i+1)] * x[k];
  37.         v[i] = -sum;
  38.         if (i != NPOLE - 1)  {
  39.             limit = NPOLE - (i+1);
  40.             for (l=0; l < limit ;++l)  {
  41.                 sum += x[NPOLE-(i+1)-(l+1)]* x[NPOLE-(l+1)] - x[FRAME-(i+1)-(l+1)]* x[FRAME-(l+1)];
  42.                 a[(i+1)+l][l] = a[l][(i+1)+l] = sum;
  43.             }
  44.         }
  45.     }
  46.     sum = (double) 0.0;
  47.     for (k=NPOLE; k < FRAME ;++k)
  48.         sum += pow(x[k], (double) 2.0);
  49.     sumy = sumx = sum;
  50.     for (l=0; l < NPOLE ;++l)  {
  51.         sum += pow(x[NPOLE-(l+1)], (double) 2.0) - pow(x[FRAME-(l+1)], (double) 2.0);
  52.         a[l][l] = sum;
  53.     }
  54.     gauss(a, v, b);
  55. /*    filtn(x, y, b);   */
  56.     for (i=0; i < NPOLE ;++i)
  57.         sumy = sumy - b[i]*v[i];
  58.     *rms1 = sqrt(sumx/(double) (FRAME - k1 + 1) );
  59.     *rms2 = sqrt(sumy/(double) (FRAME - k1 + 1) );
  60.     *errn = pow(((*rms2)/(*rms1)), (double) 2.0);
  61.     for (bp=b; bp-b < NPOLE ;++bp,++c)
  62.         *c = *bp;
  63.     return(0);
  64. }
  65.  
  66. gauss(aold, bold, b)
  67. double aold[POLEMAX][POLEMAX], *bold, *b;
  68. {
  69.     double amax, dum, pivot;
  70.     double c[POLEMAX], a[POLEMAX][POLEMAX];
  71.     int i, j, k, l, istar, ii, lp;
  72.  
  73.     istar = 0;
  74.     /* aold and bold untouched by this subroutine */
  75.     for (i=0; i < NPOLE ;++i)  {
  76.         c[i] = bold[i];
  77.         for (j=0; j < NPOLE ;++j)
  78.             a[i][j] = aold[i][j];
  79.     }
  80.     /* eliminate i-th unknown */
  81.     for (i=0; i < NPOLE - 1 ;++i)  {        /* find largest pivot */
  82.         amax = 0.0;
  83.         for (ii=i; ii < NPOLE ;++ii)  {
  84.             if (fabs(a[ii][i]) >= amax)  {
  85.                 istar = ii;
  86.                 amax = fabs(a[ii][i]);
  87.             }
  88.         }
  89.         if (amax < 1e-20)  {
  90.            /*    fprintf(stderr, "gauss: ill-conditioned.\n");*/
  91.             return(0);  /* ill-conditioned */
  92.         }
  93.         for (j=0; j < NPOLE ;++j)  {    /* switch rows */
  94.             dum = a[istar][j];
  95.             a[istar][j] = a[i][j];
  96.             a[i][j] = dum;
  97.         }
  98.         dum = c[istar];
  99.         c[istar] = c[i];
  100.         c[i] = dum;
  101.         /* pivot */
  102.         for (j=i+1; j < NPOLE ;++j)  {
  103.             pivot = a[j][i] / a[i][i];
  104.             c[j] = c[j] - pivot * c[i];
  105.             for (k=0; k < NPOLE ;++k)
  106.                 a[j][k] = a[j][k] - pivot * a[i][k];
  107.         }
  108.     }       /* return if last pivot is too small */
  109.     if (fabs(a[NPOLE-1][NPOLE-1]) < 1e-20)  {
  110.     /*    fprintf(stderr, "gauss: ill-conditioned.\n");*/
  111.         return(0);      /* ill-conditioned */
  112.     }
  113.     *(b+NPOLE-1) = c[NPOLE-1] / a[NPOLE-1][NPOLE-1];
  114.     /* back substitute */
  115.     for (k=0; k<NPOLE-1; ++k)  {
  116.         l = NPOLE-1 -(k+1);
  117.         *(b+l) = c[l];
  118.         lp = l + 1;
  119.         for (j = lp; j<NPOLE; ++j)
  120.             *(b+l) += -a[l][j] * *(b+j);
  121.         *(b+l) /= a[l][l];
  122.     }
  123.     return(1);
  124. }
  125.  
  126.  
  127. filtn(x, y, b)
  128. double x[], b[], y[];
  129. {
  130.     double sum;
  131.     int i, j;
  132.     double *yp;
  133.  
  134.     for (yp=y; yp-y < NPOLE ;++yp)
  135.         *yp = (double) 0.0;
  136.     yp = y;
  137.     for (i=NPOLE; i < FRAME; ++i)  {
  138.         sum = x[i];
  139.         for (j=0; j < NPOLE ;++j)  {
  140.             sum = sum + b[j] * x[i-(j+1)];
  141.         }
  142.         *(yp+i) = sum;
  143.     }
  144.     return(0);
  145. }
  146.